Thermo-modulational interband susceptibility and ultrafast temporal dynamics in 

nonlinear gold-based plasmonic devices 
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Starting from first principles, we theoretically model the nonlinear temporal dynamics of gold- 
based plasmonic devices resulting from the heating of their metallic components. At optical fre- 
quencies, the gold susceptibility is determined by the interband transitions around the X, L points 
in the first Brillouin zone and thermo-modulational effects ensue from Fermi smearing of the elec- 
tronic energy distribution in the conduction band. As a consequence of light-induced heating of 
the conduction electrons, the optical susceptibility becomes nonlinear. In this paper we describe, 
for the first time to our knowledge, the effects of the thermo-modulational nonlinearity of gold on 
the propagation of surface plasmon polaritons guided on gold nanowires. We introduce a novel 
nonlinear Schrodinger-like equation to describe pulse propagation in such nanowires, and we predict 
the appearance an intense spectral red-shift caused by the delayed thermal response. 
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I. INTRODUCTION AND MOTIVATIONS 

The design and development of subwavelength pho- 
tonic devices with metallic components has become a 
subject of intense research in the last decade. This trend 
is justified by the need for compact high-performance op- 
tical devices and is mainly driven by the enormous tech- 
nological improvement in nano-fabrication techniques. 

These state-of-the-art manufacturing tools for metallic 
nano-circuits have made it possible to design and engi- 
neer the effective optical properties of artificial materi- 
als, commonly known as metamaterials [U, H[. I n these 
synthetic materials, the propagation of light is strongly 
influenced by the geometric properties of the embedded 
metallic nano-circuits. In particular, the metallic nano- 
structures can be tailored in such a wa y th at the effective 
refractive index becomes negative |3l-[lll|. Negative in- 
dex materials (NIMs) are potentially important in super- 
lensing [l2l - [l4| and cloaking applications [T5MT7I ■ Novel 
physical mechanisms occur in anisotropic metamaterials, 
which under some circumstances can exhibit hyperbolic 
dispersion [l8l - [22l |. In the nonlinear regime, metamate- 
rials have been studied for harmonic generation [23M27j j . 
soliton pr opagati on [2814321 and optical modulation and 
switching [33-37]. 

The unusual properties of metamaterials arise from the 
fact that the metallic nano-circuits are much smaller than 
the wavelength of light, resulting in a space-averaged 
macroscopic dielectric response. Conversely, in the case 
where the optical wavelength is comparable with the di- 
mensions of the metallic sub-structures, the light feels the 
geometric details and plasmon polariton modes are ex- 
cited. Surface plasmon polaritons (SPPs) are electromag- 
netic waves propagating on metallic surfaces [Hj]. They 
constitute the best candidates for manipulating light 
on the nanoscale and for the development of subwave- 
length all-optical devices [39l - l42| . In particular, plas- 
monic waveguides have important applications as optical 



interconnects in highly-integrated optoelectronic devices 
[43l | . Other relevant applications of SPPs are found in 
medicine Q, sensing p5H47| and nano-lasers [48l - [54j |. 
The nonlinear properties of SPPs can be used for second 
harmonic generation (SHG) [55[ , active control [H, H3] 
and nanofocusing [H3, HI] . Nonlinear self-action can be 
exploited for manipulating transverse spatial diffraction 
by self-focusing [601 ] and for the formation of plasmon- 
solitons |6lM63| . Fundamental studies of metamateri- 
als and SPPs are closely related. Indeed, relevant phe- 
nomena occurring in metamaterials are observed also in 
plasmonics, e. g. n egative refraction [64],[65|, anomalous 
diffraction [66l |67| and electromagnetic cloaking [68W701 ] . 
In both fields, the innovative step is the use of nanostruc- 
tured metals for manipulating light. 

In most of the nonlinear studies reported above the op- 
tical response of metals is assumed to be linear, while the 
nonlinearity originates from the dielectric medium; how- 
ever, experimentalists know well that the Kerr nonlin- 
earity of metals can be enormous. Experimental results 
indicate strong third-order nonlinear susceptibilities that 
vary by several orders of magnitude, with values of X3 1 
that vary between 1CT 14 and 10 _18 m 2 /^ 2 0iZl] and 
that are much bigger than the third order susceptibil- 
ity of bulk silica (xf 4 ~ lCT 22 m 2 /^ 2 ). Recently, the 
nonlocal ponderomotive nonlinearity for a plasma of free 
electrons has been proposed as a possible model for the 
interpretation of experimental results [77, 78]. The pre- 
dicted value for the ponderomotive third-order suscepti- 
bility at optical frequencies (x3 ~ 10~ 20 m 2 /V 2 ) is how- 
ever insufficient to explain the experimental findings. In 
addition, the spectral dependence of the ponderomotive 
nonlinearity (x3 oc 1/uj 4 ) does not fit with the enormous 
spectral variation (by several orders of magnitude) ob- 
served in the measurements, suggesting that the basic 
nonlinear mechanism for metals is resonant. Theoretical 
and experimental confirmation of this hypothesis is to 
be found in the results of Rosei, Guerrisi et al. on the 
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thermo-modulational reflection spectra of thin films of 
noble metals [79T483T ] . In their work, the authors theoret- 
ically predict and experimentally observe a strong mod- 
ulation in the reflection spectrum due to light-induced 
heating. They demonstrate that the temperature change 
smears out the energy distribution of the conduction elec- 
trons, affecting the resonant interband absorption and 
hence the dielectric susceptibility. This process is intrin- 
sically nonlinear, since the temperature change modulat- 
ing the dielectric response depends on the optical power. 
Subsequent pump-probe experiments in thin films (8^ - 
l8r| and nanoparticles [13, IsH have confirmed the initial 
results of Rosei and Guerrisi. Theoretical and experimen- 
tal investigations on the temporal dynamics of the system 
clearly indicate that the nonlinear response of metals is 
characterized by a delayed mechanism [84l 89] , as is typ- 
ical for thermal nonlinearities (90j . 

Very recently, a complete analysis of the nonlinear op- 
tical response of noble metals, leading to the first theoret- 
ical derivation of a consistent model for the third-order 
nonlinear susceptibility of gold, was reported [9l|. Al- 
though experiments in thin films have been satisfactorily 
explained [91], a theoretical description of the thermo- 
modulational interband nonlinearity for ultrashort opti- 
cal pulses propagating in plasmonic waveguides is still 
missing. 

In this manuscript we derive the thermo-modulational 
nonlinear susceptibility reported in (9l| , starting from the 
band structure of gold, and describe its effect on SPPs 
propagating in a gold nanowire surrounded by silica glass. 
The paper is organized as follows. In section I we describe 
the optical properties of gold, the interband transitions, 
their effect on the dielectric susceptibility and its temper- 
ature dependence. In section II we model the temporal 
dynamics of the electrons through the two-temperature 
model (TTM), deriving the characteristic temporal re- 
sponse function. Finally, in section III we model the 
propagation of SPPs along a gold nanowire by introduc- 
ing a novel nonlinear Schrodinger-like equation and pre- 
dicting for the first time to our knowledge a strong red- 
shift caused by the thermo-modulational nonlinearity of 
gold. 



II. OPTICAL PROPERTIES OF GOLD 

The nonresonant optical properties of metals can be 
described through the free-electron model, where elec- 
trons are considered as free charges moving in response 
to an optical field Re[i?oe~ 1 "*] oscillating at angular fre- 
quency uj. In this model, the dielectric response of the 
plasma can be derived directly from the non-relativistic 
single-particle equation of motion (9^ : 
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where lo p = y / ne 2 /eom e is the plasma frequency, eo is 
the vacuum permittivity, n is the electron number den- 



FIG. 1: (a) Real and (b) imaginary parts of the dielectric 
constant of gold. The full lines represent the free-electron 
prediction ei ntr a, while the dashed lines correspond to a fit to 
the experimental data in Ref. ■ The open circles re pre sent 
the experimental data points of Johnson and Christy [941 ]. 



sity, e, m e are the electron charge and mass and 7 is a 
characteristic frequency accounting for electron-electron 
collisions. This model is justified by the fact that for 
metals the Fermi energy lies within the conduction band 
and many accessible states exist for the electrons. From a 
quantum perspective, free-electron motion only accounts 
for intraband transitions. 

For wavelengths in the far-infrared, the free-electron 
model provides very good quantitative a gree ment with 
experimental data for all noble metals 93]. For the 
special case of silver, the nonresonant model also works 
well at optical frequencies. In contrast, gold and copper 
susceptibilities have properties that are more involved 
and the nonresonant model is not appropriate at opti- 
cal frequencies or in the near infrared. Indeed, inter- 
band transitions between the d-band and the conduction 
band become more important and cannot be neglected 
if one wants to model accurately the optical response 
of such metals [Uj]. The presence of interband tran- 
sitions enriches the variety of physical processes occur- 
ring in metals and lies behind the strong temperature 
dependence of the dielectric susceptibility. In what fol- 
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lows, we focus on the particular case of gold and its di- 
electric properties. In Figs. QIa,b) the real and imagi- 
nary parts of the dielectric constant of gold are plotted. 
Full lines represent the free-electron prediction ej„ tra , 
while dashed lines correspond to a fit to the experimen- 
tal data in Ref. [93|. The open circles represent the 
experimental data points of Johnson and Christy [94j . 
For the free-electron calculations we used the parameters 
Lj p = 1.1515 x 10 16 rad/sec, 7 = 8.9890 x lO^sec" 1 , ob- 
tained by fitting e"„ tra (w) to the experimental data for 
long wavelengths in the far-infrared. Note that at optical 
frequencies the measured dielectric susceptibility devi- 
ates significantly from the predictions of the free-electron 
model as a consequence of two intense absorption peaks 
at A = 300, 410 nm. Hence the actual dielectric constant 
of gold can be expressed as the sum 



t-intra Winter 

where ti n tra{w) is given by Eq. ([1]). 



(2) 



A. Interband transitions 

In what follows, we review the perturbative theory for 
the interband transitions of electrons from the d-band to 
the conduction band, describing the effective interband 
contribution to the dielectric constant of gold ei n t er (o;). 
This contribution is strongly dependent on the electron 
temperature T e . We start our analysis from the deter- 
mination of the transition probability of a single electron 
from the valence to the conduction band, calculating the 
temperature-dependent interband absorption rate as a 
function of the optical frequency. In such a derivation, 
only direct transitions are accounted for and umklapp 
processes are neglected [92| . 

Gold is characterized by a face centered cubic (f.c.c.) 
lattice structure, sketched in Fig. [5^, where the Wigner- 
Seitz primitive cell is a rhombic dodecahedron and the 
lattice constant is a — 4.08A (92J. The reciprocal lat- 
tice, depicted in Fig. [2b, is body centered cubic (b.c.c.) 
where the Wigner-Seitz primitive cell is a truncated octa- 
hedron. At optical frequencies, the interband absorption 
is resonant around the points X and L in reciprocal space 
[83|, which correspond respectively to the centers of the 
square and hexagonal facets of the truncated octahedron 
(see Fig. [5b). Notably, such points are highly symmetric 
and around them the lattice vector k can be expressed as 
the sum k = kx + k//, where k± lies on the square (X) 
or hexagonal (L) facets and k// is perpendicular to them 
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FIG. 2: (a) Real and (b) reciprocal lattice structures of gold. 



(see Fig. [2b). In addition, the Fermi surface has cylindri- 
cal symmetry around the X, L points [92] and the valence 
and conduction bands can be approximated in (kj_,k//) 
space by elliptic and hyperbolic paraboloids [H[: 



E v (k) = E 0v 
E c (k) = E 0c , 
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TABLE I: Band parameters for the X and L transitions used 
in our calculations (m e is the electron mass) [83L l95j| . 



where v,c indicate the valence and conduction bands. 
The values of the constants E Qv , E 0c , m v ±, m c ±, m v iu 
m c // for the X and L transitions that we use in our calcu- 
lations are listed in Table [U Note that, as a consequence 
of the cylindrical symmetry of the Fermi surface, the con- 
duction and valence bands around the X, L points do not 
depend on the direction of k±, but only on its modulus 
k±. The Fermi level Ep lies in the conduction band E c (k) 
and for the sake of simplicity (and without loss of gen- 
erality) we perform a constant shift of all energies and 
assume that Ep = 0. 

The valence and conduction bands E v (k±,k//), 
E c (k±, k//) are plotted as functions of the moduli kj_,k// 
in Figs. [Sk,b for the X and L points. In this figure, the 
quantities k± 7 k// are measured in terms of fcr = V3n/a, 
which is the distance between the X, L points and the 
center of the Brillouin zone (r) (92[. The upper surfaces 
correspond to the conduction bands while the lower sur- 
faces correspond to the valence bands. Note that the 
paraboloid approximation made in Eqs. (|3I4[) is accurate 
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FIG. 3: Valence and conduction bands around the (a) L and 
(b) X points as functions of the lattice moduli kx,k/i. The 
upper surfaces represent the conduction bands and the lower 
surfaces the valence bands. The quantities k±,k// are normal- 
ized to fcr = y/Siv/a, the distance between the X, L points and 
the center of the Brillouin zone (r). 



only if fcx » k// <C &r ■ Note also that every point in the fic- 
titious (k±, k//) space corresponds to a circle of radius k± 
at a distance k // from either the X or the L points in re- 
ciprocal fc-space. The quantum states of electrons in the 
valence and conduction bands are Bloch wavefunctions 



ip v<c = Q 1/2 mt: (r)exp(ifc-r), 



(5) 



where Q is the primitive cell volume and k is the Brillouin 
wavevector. Fermi's golden rule provides the probability 
per unit time of transitions from valence to conduction 
band M: 



E B (k) - E t (k) - hu , (6) 



Qhm 2 uj 2 

where Eq is the electric field amplitude and 



Pv,- 



(7) 



In Eq. ([6]) umklapp processes have been neglected and 
the Dirac delta-function ensures conservation of energy 
for the direct transitions (k is conserved). To obtain the 
transition rate from the initial to the final band one needs 
to sum over all available states, labeled by the k vector. 
The sum over k can be replaced by integration if the 
density of states D(k) = 2£1/(2tt) 3 is introduced. Thus, 
the transition rate per unit volume is 



W c , v = 



(2tt) 3 



d 3 kR c , v f{E v ,T e )[l - f{E c ,T e )}, (8) 



BZ 



where T e is the electronic temperature, f(E,T e ) is the 
Fermi-Dirac occupation number and the integration is 
taken over the volume of the first Brillouin zone. Note 
that the Pauli exclusion principle for every v — > c di- 
rect transition is carefully considered in the expression 
for W C}V . Indeed, the factor f[E v (k),T e ] accounts for the 
probability that the fc-state in the valence band is occu- 
pied, while the factor 1 — f[E c (k),T e ] accounts for the 
probability that the fc-state in the conduction band is 
empty. The absorbed power per unit volume Pa can be 
directly related to the imaginary part of the dielectric 
constant 19711: 



Pa = W c>v huj = -e oje^ nter \E \ 2 . 



(9) 



Hence, the imaginary part of the dielectric constant 
due to the interband transitions is explicitly given by 



-(w,T e ) 
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where we have approximated the matrix element p CtV to 
be independent of k, which is true only in the limit \k\ <!C 
kr- J c ,v{w,T e ) provides the number of available direct 
v — > c transitions responsible for interband absorption. 
For this reason, this quantity is usually named the joint 
density of states (JDOS): 



Jc,v(u, T e ) 



(27T) 



/ d 3 k \d \E c (k) - E v {k) - hoj 
Jbz L 



:f[E v (k),T e }(l-f[E c (k),T e })] 



(11) 



The real part of the dielectric constant can then be ob- 
tained directly from Eq. (|10p using the Kramers-Kronig 
relation 
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where V represents the principal value of the integral. 



B. Calculation of the JDOS and of its 
thermo-derivative 

An exact analytical calculation of the JDOS, given by 
Eq. (jlll) . is not possible to the best of our knowledge, 
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FIG. 4: (a) Contour-plots of the CECBS hyperbolae in the 
(k±, k/i) plane. Blue, cyan and red curves correspond to E — 
—2, 0, +2eV, while full and dashed curves refer to the X and 
L transitions. The full and dashed black lines correspond to 
the asymptotes of the X , L hyperbolae, (b) Contour-plots of 
the CEDS ellipses in the (k±,k//) plane. Blue, cyan and red 
curves correspond to %uj — E = 2, 3, 4 eV and full and dashed 
curves refer to the X and L transitions. 



so in this section we calculate the JDOS and its thermo- 
derivative numerically. In turn, we use the JDOSs of the 
X, L transitions to calculate the total interband contri- 
bution to the dielectric constant and its dependence on 
the electronic temperature T e . 

The calculation of the JDOS can be greatly simplified 
by introducing another quantity, the energy distribution 
of the JDOS (ED JDOS) [H, defined as 

D(E,w) = [ d 3 k6(E c -E v -hu J )5(E c -E) = 

(2n)° J 



dl 



(2tt) 3 J a \V % E v xV t E c y 



(13) 



EDJDOS as the integral 

Jcv(u,T) = [ dE{D{E,u)x (14) 

JE„„ 

xf{E-hu,T e )[l-f(E,T e )]}. 

At every constant value of conduction energy E, the 
CECBS is a hyperbola in the plane of the moduli (k±, k//) 

and correspondingly a hyperboloid in fc-space. Some 
contour-plots of the CEBCS hyperbolae for the X (full 
curves) and L (dashed curves) transitions are depicted 
in Fig. [4^. Blue, cyan and red curves correspond to the 
energy values E = —2, 0, +2eV, indicated with arrows. 
The full and dashed black lines represent the asymptotes 
of the X, L hyperbolae. Note that the cyan curve repre- 
sents the Fermi level and that the concavity of the hy- 
perbolae depends on the sign of E — Eq c . For every fixed 
E, ui values such that E — Tllj < Eq v , the CEDS is an 
ellipse in the plane of the moduli (k±,k//) and corre- 
spondingly an ellipsoidal cap in fc-space. Some contour- 
plots of the CEDS ellipses for the X (full curves) and L 
(dashed curves) transitions are depicted in Fig. QJ). Blue, 
cyan and red curves correspond to the energy differences 
fiuj — E = 2, 3, 4eV, also indicated with arrows. The blue 
dashed line is absent since the condition for existence of 
the L ellipse is not fulfilled at fiuj — E = 2eV. 

In Eq. (|13|) , the volume integral in the reciprocal space 
has been reduced to a circuit integral over the closed in- 
tegration path A as a consequence of the Dirac delta- 
functions S (E c — E v — Tiuj) , (5 (E c — E). Hence, the in- 
tegration path A is a circle of radius k± displaced at a 
distance k// from the X, L points in the reciprocal space, 
resulting from the intersection of the CECBS (an hyper- 
boloid) and the CEDS (an ellipsoidal cap) . Such a circle 
in /c-space corresponds to the point (k±, k//) in the space 
of the moduli. The circuit integral in Eq. (fT3")) can be 
solved straightforwardly, leading to the following expres- 
sion for the EDJDOS: 



D(E,u) 



g£- 3 / 2 9[V{uj)-E] 
87r 2 07 y/V{uj) - E ' 



(15) 



where g is a degeneracy number related to the number of 
X,L points in the first Brillouin zone (g — 6, 8 for X, L 
transitions), Eq = Eq c — Eq v , 9(x) is the Heaviside step 
function and 



V = 

vh = 

£-3/2 = 



1 + m c j_/m v ^, 
^-'(Sw + ^oc- E G ), 
2m c ± / 2m v //m v ±m c // 



(16) 



m v //m c ± + m v j_m c /f 



where the line-integral is taken over the closed path A, 
given by the intersection of the constant energy of the 
conduction band surface (CECBS): E c (k) — E, with the 
the constant energy difference surface (CEDS): Ulu — E + 
E v (k) = 0. The JDOS can be expressed in terms of the 



Note that £ has the physical dimension of energy (eV), 
so that D(E, (J) has the physical dimension of the inverse 
of energy squared (eV~ 2 ). The EDJDOS D(E) of gold is 
plotted as a function of the energy level of the conduction 
band E in Fig. [S] for different values of Tllo. Note that 
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FIG. 5: EDJDOS of gold as a function of energy E in the 
conduction band for X (full curves) and L (dashed curves) 
interband transitions. Blue, cyan, green and red curves cor- 
respond to Hcj — 1, 2, 3, 4 eV. 



the EDJDOS is singular at the points E = V(uj), known 
in solid state physics as Van Hove singularities |92| . 

As can be understood from Figs. H^,b, the neck of 
the CECBS hyperbola grows with energy E while for ev- 
ery fixed frequency u the vertical semi-axis of the CEDS 
ellipse shrinks as E increases. Thus, there exists a maxi- 
mum energy E max where the CEDS ellipse is either tan- 
gential to the CECBS hyperbola on its neck or the semi- 
axes of the CEDS ellipse vanish: 

E max (oj) = {<r- kHuj)6(E g - huj) + V{u>)Q{hu - E G ), 

(17) 

where k = m v ///(m c // — m v //) and a = Eq c + kEq- 
Conversely, as a consequence of the paraboloid approx- 
imation for the conduction band, the lower integration 
boundary E m i n (u>) remains arbitrary. We set the mini- 
mum energy to 



E,, 



h 2 k 2 



where ki = fcr/5 and 

m v //m c //(m v ± + m c ±) 



M 



m v ±m c // + m v //m c ± 



(18) 



(19) 



The choice of the lower integration extremum is critical. 
Indeed, by choosing a small value for fc; one neglects the 
dispersion of the valence band; on the other hand, by 
choosing a large value of hi, the constant matrix element 
and parabolic band approximations cease to be valid. 

The problem of calculating the JDOS for the X, L in- 
terband transitions of gold is then reduced to the cal- 
culation of the integral in Eq. (|14p . which unfortunately 
has no simple analytical solution. By labeling the JDOSs 
of the X, L transitions with J* v , J^ v , the total interband 
correction to the imaginary part of the dielectric constant 



FIG. 6: Results of a numerical calculation of the imaginary 
part of the dielectric constant of gold 4 as a function of 
the wavelength A. The open circles are the experimental 
data points of Johnson and Christy [94l ]. The red dashed 
curve corresponds to a fit to the experimental data in Ref. 
[93| , while the black dashed curve is the intraband contribu- 
tion given by Eq. JTJ). The full blue, cyan and red curves 
correspond to the sum of the interband and intraband con- 
tributions e[Ji = € intra + e 'i'nter f° r electronic temperatures 
T e = 300, 700, 1000 °K. 



of gold is given by the sum 

7T6 2 

4t fi r("> T e) = 3£om 2 w 2 [ 



2 jX 



\P, 



2 jL 



(20) 



The real part of the interband dielectric suscepti- 
bility e' inter (uj, T e ) can be calculated from Eq. (H 



We have computed J^ V ,J^ V and e" nter (ui,T e ) numeri- 
cally. The results of these calculations are plotted in 
Fig. [S] For the dipole matrix elements |j?;fj 2 , \Pc, v \ 2 
we have used the values in Refs. [H, [9!| {9l\Pc,v\ 2 = 
1.6015 x 10- 47 Jxkg,g x \PcJ 2 = 0.321 x g L \PcJ 2 , where 
9l = 8,gx = 6). In Fig. [51 the numerical calculation of 
e m = e intra + e inter 1S plotted at several electronic tem- 
peratures T e = 300, 700, 1000 °K (blue, cyan and red full 
curves) . The open circles represent the experimental data 
of Johnson and Christy [94j , the red dashed line is a fit to 
the experimental data in Ref. [93| and the black dashed 
line represents the intraband contribution eintra- Note 
that the numerical results fit quite well to the experimen- 
tal data of Johnson and Christy for 400nm < A < 1/xm. 
In particular, a very good fit to the measurements is ob- 
tained for T e = 700 °K. For A < 410nm, other interband 
transitions become important and the contribution of the 
X, L points is not sufficient to explain the experimental 
measurements. 

If the conduction electrons are taken out of equilib- 
rium by light-induced heating, the interband absorption 
is affected by the so called Fermi smearing effect [831 ] . In- 
creasing temperature broadens the electron distribution 
around the Fermi energy, modifying the effective opti- 
cal properties of the metal. In order to understand the 
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FIG. 7: Thermo-modulation of the dielectric constant of gold, 
(a) Real and (b) imaginary corrections to the dielectric con- 
stant Ae m = d Te e inter (uj,T )(T e - T ) (T = 300 °K). Blue, 
cyan and red curves correspond to the electronic temperatures 
T e = 320, 350, 400 °K. 
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Tj x 10 7 


-1.6969 


-2.9413 


+5.0681 


-1.0016 


+0.4045 




0.3982 


0.3541 


0.3140 


0.2587 


0.2238 


Jj/Up 


0.0217 


0.0216 


0.0173 


0.0217 


0.0130 



TABLE II: Fit parameters for the thermo-derivative 
drJLteriu), given by Eq. pljl. oj p = 1.1515 x 10 16 rad/sec 
is the plasma frequency of gold, which has been calculated by 
fitting with the free-electron model the experimental data in 
the far- infrared [93j|. 



the description of light-induced self-thermo-modulation. 
Basically, as light impinges on the gold surface, the elec- 
trons in the conduction band are heated and the dielectric 
constant is modified by the amount 



Ae m (uj) = d Te e inter (u)AT e . 



(23) 



Note that this thermo-modulational process is intrinsi- 
cally nonlinear, since the increase of temperature de- 
pends on the absorbed optical power AT e (P^). The 
spectral dependence of the complex correction Ae m (u)) 
is plotted in Figs. [7^i,b for several values of temperature 
variation AT e = 20, 50, 100 °K (blue, cyan and red curves 
respectively) . Note that the spectral dependence of both 
real and imaginary parts Ae' m (u>), Ae^(tj) is non-trivial 
and they can be either positive or negative; hence the op- 
tical absorption can increase or decrease, depending on 
the wavelength A. 



III. ELECTRON TEMPORAL DYNAMICS AND 
THE TWO-TEMPERATURE MODEL 



temperature dependence of e" nter , contained in the X, L 
JDOSs, one needs to compute numerically the thermo- 
derivatives dT a J^i, L (u>,T e ). The resulting interband di- 
electric thermo-derivative can be fitted to a series of five 
Lorentzian functions: 



E 



(21) 



where uj v — 1.1515 x 10 



16 rad/sec is the plasma frequency 
of gold, calculated by fitting to the experimental data in 
the far-infrared 93] . The fit parameters are given in Ta- 
ble ITT] If a simplified Lorentzian behaviour is assumed, 
the calculation of the real part dr e e' inter through the in- 
tegration of the Kramers-Kronig relation given by Eq. 
(fT2"]l is straightforward: 



-J jj(u-Uj) 2 +7: 



(22) 



The exact calculation of the complex dielectric thermo- 
derivative dr e Winter (w) is of fundamental importance for 



As we have shown in the previous section, an optical 
beam impinging on a metal surface modifies the effective 
interband susceptibility by heating the electrons in the 
conduction band. This light-induced electron heating can 
be d escribed through the two temperature model (TTM) 
[lOOj , which takes account of the energy balance between 
the conduction electrons and the lattice. The electrons 
have a relatively small heat capacity and so thermalize 
through electron-electron collisions with a characteristic 
time of order r t h ~ 300/s. If one wishes to describe the 
temporal electron dynamics for ultrashort optical pulses 
(t w 100/s), it is necessary also to include the energy 
contribution of the non-thermalized electrons in the en- 
ergy balance. This can be calculated directly from the 
Boltzmann equation in the relaxation time approxima- 
tion [H], |89| . A phenomenological description of the elec- 
tron temporal dynamics can be obtained by separating 
the electron distributio n of ener gy into thermalized and 
non-thermalized parts [lOll . Il02j : 



d t N(t) = -( 7e 

C e d t T e (t) = C(T ; 

ddMt) = c(T e 



- 7l )N(t) + P A (t), 
T e ) + le N(t), 



(24) 



8 



where -Pa(0 is the mean absorbed power per unit vol- 
ume, T e (t),Ti(t) are the electronic and lattice tempera- 
tures, C e , Ci are the electronic and lattice heat capacities 
per unit volume and N(t) is the energy density stored 
in the non-thcrmalizcd part of the electronic distribu- 
tion. When an ultrashort optical pulse impinges on the 
metal, it is absorbed and transfers energy to the non- 
thermalized electrons. In turn, the non-thermalized elec- 
trons release energy density 7 e JV(t) to the thermalized 
electrons via electron-electron scattering and energy den- 
sity jiN(t) to the lattice via electron-phonon scattering. 
The non-thermalized electrons achieve thermal equilib- 
rium with a characteristic time delay of Tth = ile+li)^ 1 , 
where 7 e ,7; are the electron and lattice thermalization 
rates. Once heated by an ultrashort optical pulse, the 
thermalized electrons gradually release energy to the lat- 
tice via electron-phonon scattering, which is accounted 
for by the coupling coefficient C. Ultimately, for long 
times, the electrons reach thermal equilibrium with the 
lattice. The parameters used in the numerical calculation 
are given in Table IIIII Note that since the lattice heat 
capacity is much larger than the electron heat capacity, 
while the temporal variation of the electronic tempera- 
ture T e (t) is significant, the lattice temperature TJ(t) does 
not change significantly with time, i.e., T/(i) const. 

Note also that the electronic and lattice heat capacities 
C e ,Ci in principle depend on their respective tempera- 
tures T e ,Ti so that the extended TTM model is nonlin- 
ear. However, in the limit AT e (t) << T eq , it is possible 
to approximate C e , C\ as independent of their respective 
temperatures Setting d t N = and removing the 

equation for N(t) is equivalent to neglecting the ther- 
malization time T t h over which non-thermalized electrons 
release energy to the thermalized ones. This characteris- 
tic time is of order r t h ~ 300/s and can be neglected for 
long pulses. However, for pulses of duration r « 100/s 
such an approximation is not feasible and all the equa- 
tions must be retained in order to describe correctly the 
delayed nonlinearity. The TTM model can be solved 
straightforwardly in the Fourier domain, leading to the 
solution 



AT e (Aco) = T e (Acj) - Ti(Alo) = 



(25) 



T r T t h 



2L 

c, 



h T (Aw)P A (Au), 



Parameter 


Value 


Units 


Ref. 


C e 


2.1 x 10 4 




[103] 


a 


2.5 x 10 6 




[103] 


c 


2 x 10 16 


s- 1 


[103] 




2 x 10 12 


s- 1 


[84] 


H 


1 x 10 12 


-l 

s 


[84] 



TABLE III: Parameters of the two temperature model used 
in our numerical calculations and corresponding references. 
The electronic and lattice heat capacities are calculated for 



T ea = 300 °K. 
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FIG. 8: Plot of the two-temperature model (TTM) response, 
(a) Fourier transform hr(Auj) as a function of angular fre- 
quency shift Aw. The blue and red curves represent the real 
(h' T (Acj)) and imaginary (h'^(Auj)) parts, (b) TTM temporal 
response hr(t) as a function of time t. 



where Alu is the shift from the carrier angular frequency 
Wq of the ultrashort optical pulse and 



h T (Aio) = 



1 



[1 — iTthAui][l — iT r Aui] 



(26) 



In the expressions above we have used the parame- 
ter tv = C e Ci/[C{C e + Ci)] representing the relaxation 
time of the thermalized electrons with the lattice. In the 
temporal domain, the temperature variation AT e (t) = 
T e {t)-Ti{t) is given by 



1 f + °° 

AT e (t) = — / dAujAT(Auj) 

27r J-oc 



-iAtut 



(27) 



, 7e 11 



h T (t')P A {t-t')dt', 



where the temporal response function h<r{t) is 
6{t) 



h T (t) 



Tth - T r 



-t/Tth 



-t/T, 



(28) 



Note that causality is imposed by means of the Heaviside 
step function 9(t). Note also that in the limit Pa — > 
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the conduction electrons and the lattice are in equilib- 
rium and have the same temperature T e = T\ = T eq . 
Since the heat capacity of the lattice Ci is much greater 
than the electronic heat capacity C e , the temperature dif- 
ference AT e (t) = T e (t) — Ti(t) can be approximated by 
AT e (t) w T e (t) — T eq . The response functions in the fre- 
quency /it (Aw) and temporal /it(^) domains are plotted 
as functions of Aw and t in Figs. [Sk^b. In Fig. [Sh., the 
blue and red curves correspond to the real h' T (Au>) and 
imaginary /iJ.(Aw) parts. Note that, following the sign 
convention chosen for the exponential in the Fourier ex- 
pansion (e~ lAw *), a positive value of the imaginary part 
/i T (Aw) corresponds to loss, while a negative value cor- 
responds to gain. The temporal thermal response fiT(t), 
depicted in Fig. |BJd, is mainly characterized by a peak 
delayed in time by At sa 600/s. As a consequence of the 
delay, blue-shifted frequency components are suppressed, 
while red-shifted components are amplified, analogously 
to what happens i n sol id-core optical fibers as a result of 
the Raman effect jl04| . 
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FIG. 9: Thermo-modulational interband nonlinear suscepti- 
bility x>r as a function of the optical wavelength A. Blue 
and red curves correspond to the real and imaginary parts of 
(A). The region within which Xt can ^ e approximately 
considered to be a constant is also indicated. 



A. Thermo-modulational interband nonlinear 
susceptibility 

In the previous sections we have described the temper- 
ature dependence of the linear dielectric function of gold 
and the temporal dynamics of the conduction electrons 
heated by an ultrashort optical pulse. In this section we 
sum up the results obtained so far calculating the nonlin- 
ear susceptibility due to Fermi smearing of the conduc- 
tion electrons. The instantaneous power per unit volume 
absorbed by a metal is W(t) — E(t) ■ d t D{t), where 



D(t) = e Q 



+ 00 



^t-t^Eit') 



(29) 



eo is the vacuum permittivity and e m (t — t') is the tempo- 
ral dielectric response function. In the continuous wave 
(CW) monochromatic case, the electric and displacement 
fields are E(t) = E e- lUot and D(t) = e e m {uJo)Eoe~ l " ot , 
where ujq is the angular frequency. The mean ab- 
sorbed power can be calculated by averaging the in- 
stantaneous power over the fast oscillations e~ lu)ot [97j : 
Pa = (l/2)e woe^j(wo)|£^| 2 . Hence, inserting the power 
dependent temperature variation AT e into Eq. (|23[) . the 
nonlinear polarization is given by 



(30) 



X T (oJo) is the thermo-modulational interband nonlinear 
susceptibility: 



Xt fan) = ^e w e^(w )7 T (w ), 



1 



(31) 



where 



7TM = T r T th 



7e 71 | o 

c — 7) 



■M- (32) 



(3) 

The real and imaginary parts of Xt are plotted as 
functions of optical wavelength A in Fig. [9j Note that, 
as a consequence of the resonant interband transitions, 
the nonlinear susceptibility is strongly dispersive at op- 
tical frequencies and can be much greater (rs 7 orders 
of magnitude) than the Kerr susceptibility of bulk silica 
(xf* ~ 10 — 22 m 2 /V" 2 ). The strong frequency dispersion of 
gold dramatically changes its optical properties, as well 

(3) (3) 

as the signs of Re%y , Imxr • Note that, for wavelengths 
A > 750nm, the thermo-modulational nonlinear suscep- 

(3) 

tibility can be approximated as Xt ~ const. 

For ultrashort optical pulses, the calculation of the 
nonlinear dielectric polarization is more involved. In the 
slowly varying envelope approximation (SVEA) the elec- 
tric field can be expressed as E{t) = tp(t)e~' lu}ot n, where 
wo is the carrier angular frequency, n is the polarization 
unit vector and ip(t) is the envelope amplitude, which is 
slowly varying compared to the fast oscillations e _I " ot . 
The expression for the mean absorbed power Pa(£) in 
this non-monochromatic case includes also the contribu- 
tions of first-order dispersion |97j and is explicitly given 
by 



P A (t) = | 2woC(w )|^| 2 



d(ue' m 



duj 



+1 



,. d(we^) 



doj 



Wd t il>-il>d t il>') 



(33) 



In conclusion, the nonlinear polarization created by an 
ultrafast optical pulse can be expressed in terms of a 
double convolution integral 



r+oo r+oo 

P NL {t) = e / dt' / dt" x 
Jo Jo 



(34) 
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x lT {t')h T {t")P A {t -t'- t")E(t - t'). 

In this expression, 7t(0 is the interband response func- 
tion (measured in the units m 3 W~ 1 s~ 1 ) and is given by 
the inverse Fourier transform of 7t(w) (measured in the 
units m W ), which is in turn given by Eq. (|32p . 



IV. THERMO-MODULATIONAL NONLINEAR 
DYNAMICS IN PLASMONIC DEVICES 

In the previous sections we have calculated the thermo- 
modulational interband nonlinear susceptibility of gold 
starting from the basic properties of its band structure. 
In this section we study the optical propagation of surface 
plasmon polaritons (SPPs) guided along gold nanowires 
surrounded by silica glass, including the novel nonlinear 
effects originating from the heating of gold. A common 
theoretical approach to modelling optical propagation in 
optical fibres and plasmonic waveguides uses the nonlin- 
ear Schrodinger equation for the slowly varying ampli- 
tude of a guided pulse perturbatively derive d based o n 
the assumptions of low loss and nonlinearity [l04| - |lll j . 

We start the analysis from the time-dependent 
Maxwell equations for the optical electric (E) and mag- 
netic (H) fields: 



VxE(f,t) = -p d t H(r,t), (35) 
VxH(f,t) = d t D L {r,t)+d t P NL {r,t), (36) 

where r is the position vector, t is the temporal variable, 
fio is the magnetic permeability of vacuum, P/vx^i) is 
the nonlinear dielectric polarization and Dl (?, t) is the 
linear part of the displacement vector, which is given by 
the constitutive relation 



r +00 



D L (r,t) = e / dt'e L (r,t')E{r,t-t'). 



(37) 



£i(r, t) is the position-dependent temporal response func- 
tion, which can be expressed in terms of a Fourier expan- 
sion of the linear dielectric profile ei,{f, w): 



e L (f,t) 



1 

2^ 



(38) 



As a consequence of the cylindrical symmetry of the gold 
nanowire around the z-axis, the dielectric profile depends 
solely on the modulus of the position vector p = \ f\: 



£l{p, = e m (u)0(r - p) + e d (u)9(p - r), 



(39) 



where 9(x) is the Heaviside step function, lj is the an- 
gular frequency, r is the radius of the nanowire and 
e m (w),ed(w) are the linear dielectric constants of gold 
and silica. In the calculations below we use the Sellmeier 
expansion for the dielectric constant of silica ed(ui) and a 
Lorentzian fit to the experimental data for the linear di- 
electric constant of gold e m (w), as reported in Ref. [93| . 
A sketch of the plasmonic structure discussed in this sec- 
tion is depicted in Fig. [TOT 




FIG. 10: Gold nanowire of radius r surrounded by silica glass. 



A. Linear modes 

If one neglects the nonlinear polarization Pnl, Eqs. 
(|35l36p can be directly solved using the Ansatz: 



H{f,t) = e c/ 1/2 V^(p)e m * +l/3 



(40) 



where c is the speed of light in vacuum, ip is the mode 
amplitude (measured in VF 1 / 2 ), e, h are the linear guided 
mode profiles (dimensionless) , /3 is the mode propaga- 
tion constant, <f> is the angle between the vectors r, x and 
m is the azimuthal mode order. The factor J 1 / 2 is a 
constant, chosen in such a way that \ip\ 2 represents the 
total optical power carried by a linear mode with angular 
frequency u and propagation constant /3. The linear dis- 
persion relation can be directly calculated by sub- 
stituting Eqs. (|40I41I) into Eqs. (|35I36|) and by applying 
the boundary conditions for the continuity of the tan- 
gential components of the electric field and the normal 
component of the displacement vector |lllUll5 |. The 
modal profiles e, h are combinations of modified Bessel 
and Hankel functions of different orders [nHil, the 
solutions E(r,t),H(r,t) corresponding to SPPs [38[. 

Figs. Illa .b show the complex dispersion relations for 
the (a) m — and (b) m = 1 guided SPP modes; real and 
imaginary parts of the propagation constant j3 are plot- 
ted as functions of the optical wavelength A. For both the 
to — 0, 1 modes, if A > 50077,m, the real (Re/3) and imag- 
inary (Im/3) parts of the propagation constant increase 
as the optical wavelength A decreases, reaching maxima 
at A = \ sp ~ 500nm, at the surface plasmon resonance. 
For A < 500nm, the behaviour of the dispersion relation 
is more involved. In the ideal case where the metal is 
lossless, the imaginary part of the propagation constant 
vanishes, Im/3 = 0, while the real part Re/3 diverges as 
A — > A sp . 

Optical confinement depends mainly on the real part of 
the propagation constant Re/3: h igh value s of Re/3 corre- 
spond to tightly confined modes jll2l . Ill5| . On the other 
hand, the attenuation coefficient of the linear modes is 
directly related to the imaginary part of the propagation 
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FIG. 11: Linear dispersion relations /3(A) for the (a) m = 
and (b) m = 1 plasmon polariton modes. Blue, green and 
red curves correspond to the wire radii r = 50, 250, 500 nm, 
respectively. 
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FIG. 12: Contour-plots of the time-averaged Poynting vector 
for the (a) m = mode and (b) a superposition of the m = ±1 
modes on a gold nanowire with radius r = 50nm surrounded 
by silica. The optical wavelength is A = 800nm. 



constant: a = 2Im/3. Hence, if one wants to achieve tight 
SPP confinement it is impossible to avoid high losses so 
that the use of materials with large gain is required in 
practical applications [4§l-l54j]. Note that the fundamen- 
tal mode m ~ is TM polarized and that both the 
real and imaginary parts of the propagation constant 
Re/3, Im/3 increase as the radius of the gold nanowire 
decreases (see Fig. [TTk). A contour-plot of the time- 
averaged Poynting vector S z = (l/2)i • Re(E x H*) of 
the TM plasmonic mode (m = 0) of a gold nanowire with 
radius r = 50nm at optical wavelength A = 800nm is de- 
picted in Fig. H2h . Note that the electromagnetic field 
is tightly bound to the metal surface and that the power 
distribution is cylindrically symmetric. 

In contrast to the TM fundamental mode, the m = 1 
mode is hybrid polarized and less well confined. The 
dispersion relation does not depend on the sign of the 
azimuthal mode order m |l!2| so that the m — ±1 modes 
(characterized by opposite chirality) are degenerate. A 
contour-plot of the time-averaged Poynting vector S z of 
the superposition of m — ±1 guided SPP modes is shown 
in Fig. IT2"b (for the same parameters of Fig. IT2"k ) . 



Note that the power distribution of such mode is not 
azimuthally symmetric and depends on the angle <p. If 
one calculates the time-dependent Poynting vector of the 
m = ±1 modes, finds that these SPPs spiral around 
the surf ace of the gold nanowire with opposite chiral- 
ity 115j . However, the azimuthal spiralling is averaged 



out in time and as a result there is no net angular flow of 
time-averaged power for the single m = ±1 modes and 
for their superposition. In Fig. II 31 Im/3 for the m = 1 
mode is plotted as a function of r. Red, green and blue 
curves correspond to A = 700,800,900 nm. Both Re/3 
and Im/3 depend on r in a manner more complicated 
than for the fundamental TM mode (see Figs. [TTb.b lT51) . 
In particular, for fixed optical wavelength A, Im/3 is max- 
imum at a characteristic wire radius ro, decreasing sig- 
nificantly when r < tq. Hence, if the wire radius is much 
smaller than the optical wavelength r <C A, long range 
surface plasmon polaritons can be excited [3S|, Ill5j . The 
field penetration within the gold nanowire is limited for 
these modes, and hence the attenuation is significantly 
reduced. In the following nonlinear analysis we focus on 
the m = 0, 1 modes. Higher order modes are less well 
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FIG. 13: Imaginary part of the propagation constant 7m/3 as 
a function of the wire radius r for the hybrid m = 1 plasmonic 
mode. Red, green and blue curves correspond to the optical 
wavelengths A = 700, 800, 900 nm, respectively. 



confined and cut off for wavelengths greater than a par- 
ticular value X rn . 



B. Generalized nonlinear Schrodinger equation 

If the effect of loss and nonlinearity on the fast linear 
oscillations is weak, the nonlinear propagation of an op- 
tical pulse in a plasmonic waveguide can be described by 
the generalized nonlinear Schrodinger equation (GNLSE) 
for the field amplitude in the sl owly varying envelope 
approximation (SVEA) piim] . In this approach the 
Ansatz for the electromagnetic field is 

E(f,t) = I l ' 2 ij}{z,t)e{p)e im ' l>+i ^ z - iUot , (42) 
H(r,t) = e cI 1/2 tP{z,t)h(p)e m,l>+lf3oZ - lulot , (43) 

where ip(z,£) is the slowly varying envelope amplitude 
and ujq is the carrier angular frequency. With /3q we 
denote the linear propagation constant at the carrier fre- 
quency calculated by neglecting the nonlinear polariza- 
tion and the metal loss; e(p),h(p) are the corresponding 
unperturbed linear mode profiles (dimensionless) and I 
is a constant chosen so that \ip\ 2 represents the optical 
power, as in the previous section. For a gold nanowire 
surrounded by silica glass, the nonlinear polarization is 

Pnl{E) = P$l(E)9(r -p) + P^ L (E)9(p - r), (44) 

where 



P. 



Si 



\E\ 2 E + -E 2 E* 



(45) 



= 2.25 x 10 22 m 2 /V 2 is the Kerr coefficient of silica 
glass at A = 800nm and P$t{E) is given by Eq. (1341) . 
Note that in Eq. (|4"5"|) we have neglected the Raman effect 



[116j . which should in principle be retained. However, as 
we will show, the effective nonlinear coefficient of silica 
is much smaller than the effective nonlinear coefficient 
of gold so that neither Kerr nor Raman effects play any 
significant role. Hence, for the sake of simplicity, we do 
not consider the Raman term, comparing our results only 
with the Kerr term. 

Since we are mainly interested in the thermo- 
modulational interband nonlinearity, we neglect the dis- 
persive terms in the absorbed power P^t), given by Eq. 
(f3"3"|) . These terms are expected to play only a minor 
role, since they are small corrections to the carrier term. 
In what follows, we will focus on the spectral region 
(A « 800nm) where the interband response function is 
approximately constant (77^ (u>) ~ Jt(wq), see Fig. [5]) so 
that the nonlinear polarization of gold can be approxi- 
mated by 



Pnl(E) « e oX ( Z("o) I dt'h T (t')\E(t t')\ 2 E(t), 
Jo 

(46) 

where h T (t'), XaIM are given by Eqs. (|28I31I) . By 
inserting Eqs. (|42l43p into Maxwell equations and devel- 
oping a first order perturbative theory jlll| one obtains 
the GNLSE for the slowly varying amplitude ip(z,t): 

id z ip(z,t) + D(id t )i>(z,t) + r Si \ip(z,t)\ 2 ip(z,t) + 
r+00 

+T Au dt'h T {t')\i>{z,t-t')\ 2 i J {z,t)=0, (47) 
Jo 



+00 



where 



■Si 



■ Au 



(3) 
u 0X S i 



J^d0j r + °°dpp[2\e\4 



ia 2 i 



4e n c : 



ex h* 



° C (Jo" # Io + °° dppR* 

fp 2 * d4> dpp\e\ 4 



2 ' 



(3) 

U0XAu 



(lo* dfi Io + °° dppRe 



ex h* 



The linear dispersion operator D(idt) is complex, ac- 
counting as it does for the linear losses of gold. Its ac- 
tion on the envelope amplitude can be calculated in the 
Fourier domain: 



D(id t )^(z,t) 



1 

2^ 



+00 



dujD(uj)ip(z, uj)e 



where 



D(u) = P{u) 



dw 



(LU - W ). 



(48) 



(49) 



Note that flo is the real-valued carrier propagation con- 
stant of the linear unperturbed mode, /3(uj) is the com- 
plex modal wavevector calculated in the previous section 
and the prime superscript in the equation above indicates 
the real part [fi 1 = Re/3). 

The nonlinear parameters TgjjT^ are measured in 
m~ 1 W~ 1 and account also for the surface nonlinearity 
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FIG. 14: Kerr nonlinear coefScient (Tsi) of a gold nanowire 
surrounded by silica glass for the m — (full lines) and m = 1 
(dashed lines) modes. Blue and green lines correspond to the 
wire radii r — 50, 250 nm, respectively. Note that the plot is 
made in semi- logarithmic scale. 



110 L llllj , which is neglected in the averaging approach 
1041 ]. Note that while T$i is a real quantity, Tau is 
complex and accounts for the nonlinear loss of gold. The 
nonlinear parameters of silica (Tsi) and gold (Ta u ) are 
plotted as functions of the carrier wavelength Ao in Figs. 
1141151 In both figures, the full and dashed curves repre- 
sent the to = and m = 1 modes, while blue and green 
colors correspond to wire radii r = 50, 250 nm. Note 
that the real part of the gold nonlinear parameter is much 
greater than the Kerr nonlinear parameter of silica in the 
spectral region considered. Also, if r <C A, the nonlinear 
parameters of the to = 1 mode are much smaller than 
those for the m — mode since they are much less con- 
fined. In this limit, as discussed in the previous section, 
while the fundamental m — mode is tightly confined 
to the metal surface and propagates only for a few wave- 
lengths, the to = 1 mode is much less localized and can 
propagate for longer distances (long-range guided SPP 
mode). This reduction in loss is at the cost of a weaker 
effective nonlinearity. 

The formulation of the propagation equation Eq. (|47|) 
constitutes the main result of this paper. We have nu- 
merically solve d Eq . (|4"7|) using the fast Fourier split- 
step algorithm 104j . As we have already shown, when 
Aq ~ 800nm the Kerr nonlinearity of silica does not play 



a significant role and can be neglected. For the numerical 
simulations we have considered a hyperbolic secant input 
pulse: ip(0,t) — Pi n sech(t /to) , with to — 106/s. Pj„ 
is the instantaneous pulse power, which can be directly 
calculated from the average power of the laser source: 
Pin = C e ffP av /(2u rep to), where C e // is the launch effi- 
ciency of the laser beam into the gold nanowire, v rep is 
the repetition rate and 2tg is the pulse duration. The in- 
stantaneous power is kept well below the damage thresh- 
old power of gold, above which m eltin g, ablation and 
vaporization occur (P dam « 10 6 W) [Tl7| . In Fig. HH the 
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FIG. 15: (a) Real and (b) imaginary parts of the thermo- 
modulational interband nonlinear coefficient (T^„) of a gold 
nanowire surrounded by silica glass for m — (full curves) 
and m = 1 (dashed curves) modes. Blue and green curves 
correspond to the wire radii r = 50, 250 nm, respectively. 
Note that the plots are made in semi- logarithmic scale. 



numerical propagation along a gold nanowire with radius 
r = 50nm surrounded by silica glass is depicted for (a) 
to = 0, P in = lxl0 4 VFand (b) m = 1, P in = 5.3xl0 5 lF. 
In this contour-plot the modulus of the Fourier transform 
of the optical amplitude (|?/ ; (^,w)|) is shown. The m = 
TM mode is highly nonlinear and significant nonlinear 
dynamics can be observed even for relatively small opti- 
cal power. However, such a high nonlinearity is paid for 
by high loss, limiting the effective propagation length to 
L rj 2/iTO (see Fig. ITBk). The hybrid polarized m = 1 
mode supports long-range guided SPP modes so that 
both the nonlinear and loss coefficients are significantly 
smaller. In this case, in order to observe a strong nonlin- 
ear dynamics, it is necessary to use a considerably higher 
optical power. For the to = 0, 1 modes a signature red- 
shift indicates the presence of a thermo-modulational in- 
terband nonlinearity, analogously with the Raman effect 

ma. 

This red-shift is the natural consequence of the 
intrinsic delayed mechanism governing the thermo- 
modulational interband nonlinear susceptibility of gold. 
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FIG. 16: Nonlinear propagation of an optical pulse along 
a gold nanowire with radius r — 50nm surrounded by sil- 
ica glass for: (a) m = and an instantaneous input power 
Pin = 1 x 10 4 VK; (b) m = 1 and an instantaneous input power 
Pin = 5.3 x lO 5 !! 7 . In both figures the input pulse is a hyper- 
bolic secant t/>(0, t) = V Pi„sech(t /to) , with to = 106/s. The 
contour-plot displays the modulus of the Fourier transform of 
the optical amplitude: |i/>(z,u)|. 



FIG. 17: Transmission spectrum (T = ln\tp(L, w)/^om|, 
where ipoM = max[ip(0, oj)]) of the m = 1 long-range mode, 
(a) T is calculated at a propagation length of L = 100/im, for 
a fixed wire radius of r — 50nm and for different input optical 
powers: P m = 5.3 x 10 4 W (bhw curve), P ln = 2.7 x 10 5 W 
(green curve) and Pin = 5.3 x 10 J W (red curve), (b) Trans- 
mission spectrum (T) for a propagation length of L = 20/im, 
fixed input power P; n = 5.3 x 10 4 W and for different nanowire 
radii: r — 50nm (blue curve), r = lOOnm (green curve) and 
r = 400nm (red curve). The black dashed curves represent 
the input spectrum {ln\ip{Q, uj) /%1>om\)- 



In the time domain the frequency red-shift is accompa- 
nied by a small pulse delay of order w 1/s. We emphasize 
that neither the Kerr nor the Raman nonlinearities of sil- 
ica are large enough to produce the reported red-shift for 
the propagation lengths considered. The strong red-shift 
is accompanied by a large time-delayed nonlinear loss, as 
can be understood from Figs. [T7k.b. 

In Fig. 117a . the transmission spectrum (T = 
ln\ip(L, uj)/iPom\, where ipoM = max[ip(0, a;)] and L — 
100/ito) of the m = 1 long-range mode of a gold nanowire 
with radius r = 50nm is depicted for several input pow- 
ers: P m = 5.3 x 10 4 VF (blue curve), P in = 2.7 x 10 5 W 
(green curve) and Pi n = 5.3 x 10 5 W (red curve). The 
black dashed curve corresponds to the normalized in- 
put spectrum on a logarithmic scale (ln\ij){Q,uj)/il}n,M\)- 
Note that as the input power increases the red-shift in- 
creases and the transmission peak decreases accordingly 
as a consequence of nonlinear loss. Also, as the power 



increases, the transmission spectrum displays some weak 
oscillations, which resemble Kerr-related self-phase mod- 
ulation. Indeed, the dispersion length is much longer 
than the nonlinear length so that only the linear loss af- 
fects the nonlinear dynamics. 

In Fig. 117b . the transmission spectrum (T = 
ln\ip(L, aj)/-0oA/|, where L = 20/im) of the m = 1 
long-range mode is shown for several nanowire radii: 
r = 50nm (blue curve), r = lOOnm (green curve) and 
r = 400nm (red curve). The black dashed curve repre- 
sents the input spectrum (ln\?p(0, w)/?/>om|) and the input 
power is fixed to Pt n = 5.3 x 10 4 1F. Note that the linear 
properties of the m = 1 mode are non-trivial (see Fig. [T3"|) 
and, as a consequence, the power-dependence of the red- 
shift and the transmission peak is non-monotonic. This 
means that an optimal radius r — r Q exists, where the 
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r (nm) 

FIG. 18: Wavelength red-shift AA as a function of the wire 
radius r for the m = 1 long-range mode, fixed input power 
Pin = 5.3 x lO 4 !! 7 and the propagation length L — 20fim. 



achievable red-shift is maximum. In Fig. [T51 the red- 
shift of the transmission peak of the m = 1 long-range 
mode is plotted as a function of wire radius, for fixed 
input power Pj„ = 5.3 x 10 4 PF and for a propagation 
distance L = 20(im. The open blue circles represent the 
results of numerical simulations, while the black curve 
corresponds to an interpolation of the numerical results. 
The maximum red-shift attainable at this input power is 
AAw 7nm for a radius of r a sa llOnm. The dependence 
of the thermo-modulational interband nonlinear coeffi- 
cient (T Au) and the absorption coefficient (a — 2/3'q) on 
r strongly affects the nonlinear dynamics and as a con- 
sequence the red-shift attainable, which reaches a maxi- 
mum at r — r a . If the input power increases, the maxi- 
mum red-shift increases accordingly. 

V. SUMMARY 

In this paper we have described the thermo- 
modulational interband nonlinearity of gold starting from 
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its band structure. Electrons in the conduction band are 
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metal nonlinearity, which we have found to be several 
orders of magnitude larger than the Kerr nonlinearity of 
fused silica. We have derived, for the first time to our 
knowledge, a generalized nonlinear Schrodinger equation 
suitable for modeling the optical propagation of SPPs 
along a gold nanowire surrounded by silica glass. Solv- 
ing this equation using a fast Fourier split step algo- 
rithm, we have found that the signature of the thermo- 
modulational interband nonlinearity is a red-shift of the 
optical pulse. This red-shift results from the intrinsic 
time-delayed nature of the thermo-modulational inter- 
band nonlinearity of gold. We believe that this novel 
nonlinear effect may be important for frequency conver- 
sion in plasmonic devices. We have also provided some 
details on the expected red-shift, its dependence on the 
wire radius and the optical power necessary to observe it 
experimentally. 
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